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Abstract 

This paper pursues the study carried out by the authors in Stability and 
Hopf bifurcation in a hexagonal governor system [13], focusing on the codi- 
mension one Hopf bifurcations in the hexagonal Watt governor differential 
system. Here are studied the codimension two, three and four Hopf bi- 
furcations and the pertinent Lyapunov stability coefficients and bifurcation 
diagrams, illustrating the number, types and positions of bifurcating small 
amplitude periodic orbits, are determined. As a consequence it is found 
an open region in the parameter space where two attracting periodic orbits 
coexist with an attracting equilibrium point. 
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1 Introduction 



The centrifugal governor is a device that automatically controls the speed of 
an engine. The most important one is due to James Watt -Watt governor- 
and it can be taken as the starting point for automatic control theory. Cen- 
trifugal governor design received several important modifications as well as 
other types of governors were also developed. From MacFarlane [Sj, p. 251, 
we quote: 

"Several important advances in automatic control technology were made in 
the latter half of the 19th century. A key modification to the flyball governor 
was the introduction of a simple means of setting the desired running speed 
of the engine being controlled by balancing the centrifugal force of the flyballs 
against a spring, and using the preset spring tension to set the running speed 
of the engine" . 

In this paper the system coupling the Watt governor with a spring (resp. 
Watt governor) and the steam engine will be called simply the Watt Governor 
System with Spring (WGSS) (resp. Watt Governor System (WGS)). The 
stability analysis of the stationary states and small amplitude oscillations of 
this system will be pursued here. 

The first mathematical analysis of the stability conditions in the WGS was 
due to Maxwell [9] and, in a user friendly style likely to be better understood 
by engineers, by Vyshnegradskii [15]. A simplified version of the WGS local 
stability based on the work of Vyshnegradskii is presented by Pontryagin 

From the mathematical point of view, the oscillatory, small amplitude, 
behavior in the WGS can be associated to a periodic orbit that appears from 
a Hopf bifurcation. This was established by Hassard et al. in [5] , Al-Humadi 
and Kazarinoff in [1] and by the authors in [TT1 [T2 |. Another procedure, 
based in the method of harmonic balance, has been suggested by Denny [3] 
to detect large amplitude oscillations. 

In [TT] we characterized the surface of Hopf bifurcations in a WGS, which 
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is more general than that presented by Pontryagin [TU], Al-Humadi and 
Kazarinoff pQ and Denny [3]. 

In [12] restricting ourselves to Pontryagin's system of differential equa- 
tions for the WGS, we carried out a deeper investigation of the stability of 
the equilibrium along the critical Hopf bifurcations up to codimension 3, hap- 
pening at a unique point at which the bifurcation diagram was established. 
A conclusion derived from the diagram implied the existence of parameters 
where the WGS has an attracting periodic orbit coexisting with an attracting 
equilibrium. 

In [13] we characterized the hypersurface of Hopf bifurcations in a WGSS. 
See Theorem 14 . 1 1 and Fig. [2] for a review of the critical surface where the first 
Lyapunov coefficient vanishes. 

In the present paper we go deeper investigating the stability of the equi- 
librium along the above mentioned critical surface. To this end the second 
Lyapunov coefficient is calculated and it is established that it vanishes along 
two curves. The third Lyapunov coefficient is calculated on these curves and 
it is established that it vanishes at a unique point. The fourth Lyapunov 
coefficient is calculated at this point and found to be negative. See Theorem 
14. 2[ The pertinent bifurcation diagrams are established. See Fig. [6] and 
A conclusion derived from these diagrams, concerning the region — a solid 
"tongue" — in the space of parameters where two attracting periodic orbits 
coexist with an attracting equilibrium, is specifically commented in Section 

El 

The extensive calculations involved in Theorem 14.21 have been corrobo- 
rated with the software MATHEMATICA 5 p2] and the main steps have 
been posted in the site [16J. 

This paper is organized as follows. In Section [2] we introduce the differen- 
tial equations that model the WGSS. The stability of the equilibrium point 
of this model is analyzed and a general version of the stability condition is 
obtained and presented in the terminology of Vyshnegradskii. The Hopf bi- 
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furcations in the WGSS differential equations are studied in Sections [3] and 
HI Expressions for the second, third and fourth Lyapunov coefficients, which 
fully clarify their sign, are obtained, pushing forward the method found in 
the works of Kuznetsov [5J [7j. With this data, the bifurcation diagrams are 
established. Concluding comments, synthesizing and interpreting the results 
achieved here, are presented in Section [5j 



2 The Watt governor system with spring 
2.1 WGSS differential equations 

The WGSS studied in this paper is shown in Fig. [TJ There, ip E (0, f ) is 
the angle of deviation of the arms of the governor from its vertical axis Si, 
Q G [0, oo) is the angular velocity of the rotation of the engine flywheel D, 9 
is the angular velocity of the rotation of Si, I is the length of the arms, m is 
the mass of each ball, if is a sleeve which supports the arms and slides along 
Si, T is a set of transmission gears and V is the valve that determines the 
supply of steam to the engine. 




Figure 1: Watt centrifugal governor with a spring - steam engine system. 



The WGSS differential equations can be found as follows. For simplicity, 



we neglect the mass of the sleeve and the arms. There are four forces acting 
on the balls at all times. They are the tangential component of the gravity 

— mg simp, 

where g is the standard acceleration of gravity; the tangential component of 
the centrifugal force 

m I simp 9 2 cos<£>; 
the tangential component of the restoring force due to the spring 

—2kl{l — cosy?) simp, 

21 is the natural length of the spring and k > is the spring constant; and 
the force of friction 

—blip, 

b > is the friction coefficient. 

From the Newton's Second Law of Motion, and using the transmission 
function 9 = c Q, where c > 0, one has 



<p = h c f2 simp cos ip sirup up. (1) 

\m J ml m 

The torque acting upon the flywheel D is 

I tl = fi cos(p- F, (2) 

where I is the moment of inertia of the flywheel, F is an equivalent torque of 
the load and /i > is a proportionality constant. See [10] . p. 217, for more 
details. 

From Eq. (pQ) and ([2D the differential equations of our model are given by 

d¥ _ ^ 
dr 

dip (2k 2 2 \ . 2kl + mg . b 

—— = h c i I sm ip cos ip : sm ip ip Id) 

dr \m J ml m 

dtt 1 . 

— = - (/xcos^-F) 



where r is the time. 

The standard Watt governor differential equations in Pontryagin [10] , p. 
217, 

d 09 

dip 2 2 • 9 ■ b , , A S 
— — = c \l sin oj cos oj — - smo9 y> (A) 

dr I m 

dn l 

are obtained from Q by taking k = 0. 

Defining the following changes in the coordinates, parameters and time 

/ ml \ 1 ^ 2 ( ml \ ^ 2 

\2/a + mg ) \2kl + mg J 

2kl + mg\ 1/2 2kl 



ml J 2kl + mg 

b ( ml \ 1//2 c/i / mZ \ F 

m \2A;Z + mg ) I \2kl + mg J ' /i ' 

where 0</c<l,e>0, a > and < /3 < 1, the differential equations (j3J) 
can be written as 

x = — = y 
dt y 

y' — -jt = (z 2 + k) sinx cosx — sinx — e y (5) 

dz 

z' = — = a(cosx — 0) 
or equivalently by 

x' = /(x,C), (6) 

where 

/(x, Q = (y, (z 2 + k) sin x cos x — sin x — e y, a (cos x — (3)) , 
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and 



x= (x,y,z) e(o^)xlx [0,oo) 



C = (/3, a, e, k) e (0, 1) x (0, oo) x (0, oo) x [0, 1) . 



2.2 Stability analysis of the equilibrium point 

The WGSS differential equations (jHJ) have only one admissible equilibrium 
point 



Po = (x ,yo,z ) = [ arccos/?,0, ( ^ - k 



1/2^ 



(7) 



The Jacobian matrix of / at P has the form 

/ 1 \ 



Df (P ) 



-UJn 



-e £ 

\ -a(l -{3 2 ) 1 ' 2 j 



(8) 



where 



and 



u 



1-/3 2 



(9) 



t = 2(3 1 l 2 {l-(3 2 Y'\l-K(3) l l 2 . 

For the sake of completeness we state the following lemma whose proof 
can be found in [10], p. 58. 

Lemma 2.1 The polynomial L(X) = PqX 3 +pi\ 2 +P2X+P3, Po > 0, with real 
coefficients has all roots with negative real parts if and only if the numbers 
Pi,P2,P3 are positive and the inequality pvpi > P0P3 is satisfied. 



Theorem 2.2 If 



e > e r 



2a^ 2 {l-K(3) 1 ' 2 , 



(10) 
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then the WGSS differential equations |3P have an asymptotically stable equi- 
librium point at Pq. If 

< e < e c 

then P Q is unstable. 

Proof. The characteristic polynomial of Df (Pq) is given by p(X), where 

-p(X) = X 3 + Pl X 2 +p 2 X+p 3 , 



Pi =e, p 2 = — , Ps 



l-(3 2 2«/? 3 / 2 (l - /? 2 ) (1-k(3) 



1/2 



' "° P 

The coefficients of —p(X) are positive. Thus a necessary and sufficient con- 
dition for the asymptotic stability of the equilibrium point Pq, as provided 
by the condition for one real negative root and a pair of complex conjugate 
roots with negative real part, is given by (jlOj) . according to Lemma [2. 1[ 



In terms of the WGSS physical parameters, condition (fit)]) is equivalent 

to 

— T]>1, (11) 

m 



where 



V 



dQ 



dP 



1 (12) 



2/5 3 / 2 (i - Kpy/ 2 

is the non-uniformity of the performance of the engine which quantifies the 
change in the engine speed with respect to the load (see [TO], p. 219, for more 
details). Eq. f|T2|) can be easily written in terms of the original parameters. 

The rules formulated by Vyshnegradskii to enhance the stability follow 
directly from ( JTTT) . In particular, the interpretation of (TTTI) is that a sufficient 
amount of damping b must be present relative to the other physical parame- 
ters for the system to be stable at the desired operating speed. The condition 
( Hip is equivalent to the original condition given by Vyshnegradskii for the 
WGS (see [10], p. 219). 
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In section 0] we study the stability of P under the condition 



e = e c , (13) 

that is, on the hypersurface — the Hopf hypersurface — complementary to 
the range of validity of Theorem 12 .21 

3 Lyapunov coefficients 

The beginning of this section is a review of the method found in [Bj, pp 
177-181, and in [7] for the calculation of the first and second Lyapunov co- 
efficients. The calculation of the third Lyapunov coefficient can be found 
in [12]. The calculation of the fourth Lyapunov coefficient has not been 
found by the authors in the current literature. The extensive calculations 
and the long expressions for these coefficients have been corroborated with 
the software MATHEMATICA 5 pi]. 
Consider the differential equations 

x' = /(x,/i), (14) 

where x G K n and fi G M m are respectively vectors representing phase vari- 
ables and control parameters. Assume that / is of class C°° in IR n x IR m . 
Suppose (fl4l) has an equilibrium point x = x at ji = /i and, denoting the 
variable x — x also by x, write 

F(x) = /(x,/i ) (15) 

as 

F(x) = Ax + ~ B(x, x) + i C(x, x, x) + ^ D(x, x, x, x) + £(x, x, x, x, x) + 

— K(x, x, x, x, x, x) + L(x, x, x, x, x, x, x) + M(x, x, x, x, x, x, x, x) (16) 

H N(x, x, x, x, x, x, x, x, x) + Ofllxl 1 10 ), 

362880 v ; vm ii a 
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where A = / x (0, /xq) an d 



Ci(x,y,z) = 



5=0 



n 



d£j di k d£i 



Xj Vk zu 



e=o 



A(x,y,z,u) = V 

i>fe y=i ^' ^ ^ lKl 



5=o 



£ ' (xJ ' Z ' U,v)= ,f 

j,k,l,r,p=l 



5=0 



(17) 

(18) 
(19) 
(20) 



Ki(x,y,z,u,v,w) = ^2 



Lj(x,y,z,u, v,w,t) = ^ 

j,...,h=l 



! 90 d& di r di p d^ q 



Xj y k Z\ u r v p w q , (21) 



5=0 



d£j di k d& d£ r d£ p d£ q d£ h 



xj y k ziu r VpW q t h , (22) 



d 8 F(0 



M i( x,y,z,u,v,w,t,r) = £ 

j,...,a=l ■> 

iVj(x,y,z,u,v,w,t,r,s) = 2^ ^7 ^ 



j,...,b=l 



5=0 

y fc zi u r v p w q t h r a , (23) 

5=0 

Xj y k zi u r v p w q t h r a s b , (24) 



5=0 



for i = 1, . . . , n. 

Suppose (xo, Ho) is an equilibrium point of (JT3J) where the Jacobian matrix 
A has a pair of purely imaginary eigenvalues A2,3 = Oujq, u > 0, and admits 
no other eigenvalue with zero real part. Let T c be the generalized eigenspace 
of A corresponding to A2,3- By this is meant that it is the largest subspace 
invariant by A on which the eigenvalues are A2,3- 

Let p, q G C n be vectors such that 



Aq = iuj q, A T p = -iuj p, (p, q) = £ Pi ^ = 1, 



(25) 



i=l 
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where A T is the transposed matrix. Any vector y G T c can be represented 
as y = wq + wq, where w = (p, y) G C. The two dimensional center manifold 
can be parameterized by w,w, by means of an immersion of the form x = 
H(w,w), where H : C 2 — > W 1 has a Taylor expansion of the form 

H(w,w) = wq + wq+ V h jk w j w k + 0(\w\ 10 ), (26) 

2<j+fc<9 J 

with hj k G C n and hj k = h k j. Substituting this expression into (|T4j) we 
obtain the following differential equation 

H w w' + H & w' = F{H(w,w)), (27) 

where F is given by (|T5|) . 

The complex vectors hij are obtained solving the system of linear equa- 
tions defined by the coefficients of fl27j) . taking into account the coefficients 
of F, so that system ( 1271) . on the chart w for a central manifold, writes as 
follows 

w' = iuj w+^ G 2 \w\w\ 2 +-^ G 32 w\w\ 4 +-^ Ga3W \ w \ &+ 2^80 g 5aw\w\ 8 +0(\w\ 10 ), 

with Gj k G C. 

The first Lyapunov coefficient li is defined by 

h = g Re G 2i, (28) 

where 

G 2 i = (p, H 2 i), and H21 = C(q, q, q) + B(q, h 20 ) + 2B(q, h n ). 

The complex vector h 2 \ can be found by solving the nonsingular (n + 1)- 
dimensional system 





11 



with the condition (p, h 2 \) = 0. See Remark 3.1 of [33]. The procedure above 
can be adapted in connection with the determination of h 32 and h 43 . 
Defining H32 as 

H 32 = QB(h n , h 21 ) + B(h 20 , h 30 ) + 3B(h 2h h 20 ) + 3B(q, h 22 ) 
+2B{q, h 31 ) + 6C(g, h n , h n ) + 3C(q, h 20 , h 20 ) + 3C(q, q, h 21 ) 
+6C(q, q, h 21 ) + 6C(q, h 20 , hu) + C(q, q, h 30 ) + D(q, q, q, h 20 ) 
+QD(q, q, q, h n ) + 3D(q } q, q, h 20 ) + E(q, q, q, q, q) 

—6G 2 ih 2 i — 3G 2 \h 2 i, 

and from the coefficients of the terms w 3 w 2 in (1271) . one has a singular system 
for h 32 

(iuj I n ~ A)h 32 = H 32 - G 32 q, 
which has solution if and only if 

(p,H 32 -G 32 q) = 0. (29) 

The second Lyapunov coefficient is defined by 

Z 2 = ^ReG 32 , (30) 

where, from ([29]), G 32 = (p,H 32 ). 

The complex vector h 32 can be found solving the nonsingular (n + 1)- 
dimensional system 



A q\ 


II 


( h 32 \ 


H 


0) 




K s ) 






p 

with the condition (p, h 32 ) = 0. 
Defining 7-^43 as 

H 43 = l2B{h n ,h 32 ) + 6B(h 20 , h 32 ) + 3B{h 20 , h 41 ) 
+18B(h 21 ,h 22 ) + l2B(h 2X Ml) + 4B{h 30 , h 3 i) + B{h 30 , M 
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+4B(q, h 33 ) + 3B(q, /i 42 ) + 36C(/»n, hn, h 21 ) + 36C(h n , h 20 , h 21 
+12C(h 11 ,Ji2Q, h 3Q ) + 3C(h 20 , ^20, ^-30) + 18C(/i 2 o, ^20 , /121 
+36C(g, /in, h 22 ) + 12C(g, /i 20 , ^31) + 12C(g, /i 20 , /131 
+36C(g, h 21 ,h 21 ) + 4C(g, /i 30 , /i 3 o) + 6C(g, g, /i 32 
+12(7(5, g, M + 24C(g, /in, /i 31 ) + 18C(g, /i 20 , /i 22 
+3C(q, h 20 , h m ) + 18C(g, fo 2 l, foi) + 12C(g, ft 2 i, h 30 
+3C(q, q, h Al ) + 2AD(q, h u , hn, h u ) + 36D(g, /in, /i 20 , h 20 
+36D(g, g, hn, h 2 %) + 6L>(g, g, h 20 , h 30 ) + 18£>(g, 5, h 20 , h 21 
+4D(q, q, q, h 31 ) + l8D(q, q, q, h 22 ) + 72D{q, q, h u , h 21 
+36D(q, q, h 20 , h 2 \) + l2D(q, q, h 20 , h 30 ) + l2D(q, q, g, h 3 i 
+36D(q, /in, hn, h 20 ) + 9D(q, h 20 , h 20 ,h 20 ) + 12D(q,q,h u , h 30 
+ l8D(q, q, h 20 ,h 21 ) + D(q, q, q, h m ) + 12E(q,q,q,h u ,h 20 
+E{q, q, q, q, h 3 o) + ^E(q, q, q, q, h 21 ) + 36E(q, q, q, h u , hn 
+l8E(q, q, q, h 20 , h 20 ) + l8E(q, q, q, q, h 2l ) + 36E(q, q, q, h u ,h 20 
+4£(g, g, q, g, h 30 ) + 3E(q, q, q, h 20 , h 20 ) + 3K(q, q, q, q, q, h 20 
+12K(q, g, g, g, g, /in) + QK(q, g, g, g, q, h 20 ) + L(g, g, g, g, g, g, g 
-6(2G 32 /i 2 i + G 32 h 21 + 3G 21 h 32 + 2G 21 h 32 ), 

and from the coefficients of the terms w A w 3 in fl27j) . one has a singular system 
for /i 43 

(io; 4 - A)/i 43 = 7i 4 3 - G A3 q 
which has solution if and only if 



(P, T<43 ~ G A3 q) = 0. 

The third Lyapunov coefficient is defined by 

1 



144 



Re G 



43, 



(31) 



(32) 



where, from f[3"Tj) . G 43 = (p, 7^4 3 ). 

The complex vector /i 43 can be found solving the nonsingular (n + 1)- 
dimensional system 



iuj I n 



P 



A q\ 




( h 43 \ 


l-( 


1 


II 


V a J 
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G 4 3g 
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with the condition (p, h 43 ) = 0. 

Defining 7-^54 by the below expression 



20fl(/m, h 3 ) + lOB(h w , h 43 ) + 6B(h 20 , h 52 ) + 40B(h 21 ,h 33 ) + 30B(h 
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WB(h 22 , h 32 ) + 10B(h 3Q , h A2 ) + 4B(h 30 , h 51 ) + 40B(h 31 ,h 32 ) + 20B(h 31 , h 4 i 
5B(h 4Q , hi) + B(h 4Q ,h 5Q ) + 5B(q, h 44 ) + 4B(q, h 53 ) + 120C(h n , h n , h 32 
60C(hn,ho, h 41 ) + 360C(/i n , h 21l h 22 ) + 240C{hn,hi,hi) + 80C{h n ,h 3Q , hi 
20C(h u , ho, h 40 ) + l20C(h 20 , hn,hi) + 15C(h 20 , h 20 , hi) + 60C(h 20 , ho, h 32 
120C(/i2o, hi,hi) + l80C(h 20 , hi,hi) + lOC(h 20 , h 30 , ho) + 40C(/i 20 , ho, hi 
3C(ho, ho, ho) + l20C(ho, hi, hi) + 30C(/t 20 , hi, ho) + 60C{ho, ho, hi 
180C(h 21 ,h 21 , hi) + eOC(hi, hi, ho) + 40C(/i 30 , hi, ho) + 80C(q, h n , h 33 
30C(q, h 20 , hi) + 30C7(q, ho, h 2 ) + 120C(q, h 2U hi) + 120C(q, hi, h 3 2 
90C(q, h 22 , h 22 ) + 20C(q, h 30 , hi) + 20C{q, ho, hi) + 80C(q, h 3U hi) + 5C(q, h 40 , ho 
lOC(q, q, h 43 ) + 20C(q, q, h 43 ) + 60C{q, hn,h 42 ) + 4QC{q, h 2Q , h 33 ) + l2C{q, ho, hi 
120C(q, h 2U h 32 ) + 60C(g, hi, hi) + 40C{q, h 30 , hi) + 4C(g, ho, ho 
120C{q, h 3U h 22 ) + 20C(q, h 40 , hi) + 6C(g, q, h 52 ) + 240D(h n , h n ,hn,hi 
l2QD(hn, hn,ho, ho) + 360D(/i 20 , hi, hi, hi) + 360D(h 20 , hn,ho, hi 
60D(/i 20 , ho, hi, ho) + 90D(h 20 , h 20 , ho, hi) + 30L>(/i 20 , ho, ho, ho 
360D{q, hu,hu,h 22 ) + 240D(q, h lu ho, hi) + 720D{q, h u ,h 2U h 2 i 
80D(q, hn,ho, ho) + 240D(q, h 20 , hn,hi) + ^D(q, h 2Q , h 2Q , ho 
\80D{q, h 20 , ho,h 2 ) + 120D(q, h 20 , hi,ho) + 180D(q, h 20 , hi, hi 
15D(q, ho, ho, ho) + l80D(q, ho, ho, hi) + l20D(q, ho, ho, hi 
l20D(q, q, /in, ^32) + WD(q, q, h 20 , hi) + 60D(q, q, ho, hi) + 120D(q, q, h 2U hi 
l8QD{q, q, hi,h 2 ) + lOD(q, q, h 30 , ho) + 40£>(<7, q, ho, hi) + WD(q, q, q, h 42 
40D(q, q, q, h 33 ) + 240D{q, q, hn, h 32 ) + 120D(q, q, h 20 , hi) + 60D{q, q, ho, hi 
360£>(<7, q, h 21 ,h 22 ) + 240D(q, q, hi, hi) + 80D(q, q, h 30 , hi) + 20D(q, q, ho, ho 
30D(q, q, q, h 42 ) + 240D(q, hn,hn, h 31 ) + 60D(q, hn, ho, ho) + 360D(q, hn,h 2 i, h 21 
240D(q, hu,h 30 , hi) + 360D(q, h 20 , h n ,h 22 ) + 60D(q, h 20 , h 20 , hi 
120D(q, h 20 , ho, hi) + 360D(q, h 2Q , h 21 ,hi) + 40D(q, h 20 , h 3Q , ho 
120D(q, ho, ho, hi) + 60D(q, q, hu,h 4 i) + 60D(q, q, h 2a ,h 32 ) + 6D(q, q, h 2a , h 50 
120D(q, q, h 21 ,h 31 ) + 30D(q, q, h,ho) + Q0D(q, q, h 30 , h 22 ) + 4D(q, q, q, h 51 
120E(q, hn, hn, hn, hn) + 360£%, h 20 , hn,hi,ho) + 45S(g, h 20 , h 20 , ho, ho 
360E(q, q, hn, hi,hi) + 360£%, q, h n , ho, hi) + l20E(q, q, h 20 , hi,ho 
\8QE(q, q, h 2a , ho, hi) + 30E(q, q, ho, ho, ho) + 80E(q, q, q, hn, hi 
lOE(q, q, q, h 20 , ho) + 60E(q, q, q, ho, hi) + 40£%, q, q, h 21 ,h 3 o) + 60E(q, q, q, hi, hi 
5E(q,q,q,q, hi) + 40E(q,q,q,q, /j 32 ) + 360E(q, q, q, hn, h 2 ) + l20E(q,q,q, h 20 ,h 3 i 
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120 E(q, q, q, h 20 , h 31 ) + 360E(q, q, q, h 21 , h 2 i) + AOE(q, q, q, h 30 , h 30 ) + 60E(q, q, q, q, h 32 ) + 
720E(q, q, hu, hu,hn) + 240E(q, q, h n , h 20 , h 30 ) + 720E{q, q, h 20 , h n , h 21 ) + 
60E(q, q, h 20 , h 20 , h 30 ) + 360 E(q, q, h 20 ,h 20 , h 21 ) + 2A0E(q, q, q, hn, h 31 ) + 
180E(q, q, q, h 20 , h 22 ) + 30E(q, q, q, h 20 , h i0 ) + 180E(q, q, q, h 2 i, h 21 ) + 120 E(q, q, q, h 30 , h 21 ) + 
20E(q, q, q, q, h A1 ) + 240£'(g, h 20 , h n , h n ,h n ) + 180 E{q, h 20 , h 20 , h u , h 20 ) + 
120E(q, q, h n ,hii,h 30 ) + 360E(q, q, /i 20 , hn, h 21 ) + 90E(q, q, h 20 , h 20 , h 21 ) + 
OOE(q, q, /i 20 , h 20l h 3Q ) + 20E(q, q, q, hn,h i0 ) + -10E(q, q, q, h 20 , h 3 i) + 40£(g, q, g, /i 30 , h 21 ) + 
E{q, q, q, q, h 50 ) + 120K(q,q,q,hn,hn,h 2 o) + 30K(q,q,q,h 2Ql h 2a , h 20 ) + 
20K(q, q, q, q, hn, h 30 ) + 30K(q, q, q, q, h 20 , h 21 ) + K(q, q, q, q, q, h 40 ) + 20K(q, q, q, q, q, h 31 ) + 
2A0K(q,q,q,q,h n ,h 21 ) + AOK(q,q,q,q, h 2Q ,h 30 ) + 120K(q,q,q,q, h 2Q ,h 21 ) + 
60K(q, q, q, q, q, h 22 ) + 2A0K(q, q, q, /in, /in, /in) + 360K (q, g, q, h 20 , hn,h 20 ) + 
360K(q, q, q, q, hn, h 21 ) + 180K(q, q, q, q, h 20 , h 21 ) + 60K(q, q, q, q, h 20 , h 30 ) + 
40K(q, q, q, q, q, h 3 i) + 360K(q, q, q, h 20 ,hn, hn) + 90K(q, q, q, h 20 , h 2Q , h 20 ) + 
80K(q, q, q, q, hn,h 30 ) + 120K(q, q, q, q, h 20 , h 21 ) + 5K(q, q, q, q, q, /i 40 ) + 
60K(q, q, q, h 20 , /i 2 o, /in) + 10 K (q, q, q, q, h 20 , h 30 ) + 3L(q, q, q, q, q, h 20 , h 20 ) + 
4i(g, q, q, q, q, q, ho) + QOL(q, q, q, q, q, hn,h 20 ) + 30L(q, q, q, q, q, q, h 21 ) + 
120L{q, q, q, q, q, hn,h n ) + 60L(g, q, q, q, q, h 20 , h 20 ) + 40L(g, q, q, q, q, q, h 21 ) + 
120L(q, q, q, q, q, h 20 , h n ) + 10L(q, q, q, q, q, q, h 3Q ) + 15L(q, q, q, q, q, h 20 , h 20 ) + 
6M(q, q, q, q, q, q, q, h 20 ) + 20M(q, q, q, q, q, q, q, hu) + 10M(q, q, q, q, q, q, q, h 20 ) + 

N(q,q,q,q,q,q, q, q, q), 

and from the coefficients of the terms w 5 w 4: in fl27j) . one has a singular system 
for h 5 4 

{iu I n - A)h 54 = H 5i - G 54 q 
which has solution if and only if 

(p,H 5i -G 5 ,q) = 0. (33) 

The fourth Lyapunov coefficient is defined by 

'*=2580 ReGM - (34) 
where, from ( 1331) . G54 = (p, 7i 54 ). 

Remark 3.1 Other equivalent definitions and algorithmic procedures to write 
the expressions for the Lyapunov coefficients = 1, 2, 3, 4 ; for two dimen- 
sional systems can be found in Andronov et al. [2] and Gasull et al. j^jj, 
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among others. These procedures apply also to the three dimensional systems 
of this work, if properly restricted to the center manifold. The authors found, 
however, that the method outlined above, due to Kuznetsov FJJ [7j/, requiring 
no explicit formal evaluation of the center manifold, is better adapted to the 
needs of this work. 

A H op f point (x ,/^o) is an equilibrium point of (|14p where the Jacobian 
matrix A = /x( x o,/Wo) has a pair of purely imaginary eigenvalues A 2 ,3 = 
±ia>o, <^o > 0, and admits no other critical eigenvalues — i.e. located on the 
imaginary axis. At a Hopf point a two dimensional center manifold is well- 
defined, it is invariant under the flow generated by (TH1) and can be continued 
with arbitrary high class of differentiability to nearby parameter values. In 
fact, what is well defined is the oo-jet — or infinite Taylor series — of the 
center manifold, as well as that of its continuation, any two of them having 
contact in the arbitrary high order of their differentiability class. 

A Hopf point is called transversal if the parameter dependent complex 
eigenvalues cross the imaginary axis with non-zero derivative. In a neighbor- 
hood of a transversal Hopf point — HI point, for concision — with 1^0 the 
dynamic behavior of the system (TH|) . reduced to the family of parameter- 
dependent continuations of the center manifold, is orbitally topologically 
equivalent to the following complex normal form 

w ' = (v + iw)w + liw\w\ 2 , 

w G C, 7], u and l\ are real functions having derivatives of arbitrary high 
order, which are continuations of 0, uj q and the first Lyapunov coefficient at 
the HI point. See [6]. As li < {l\ > 0) one family of stable (unstable) 
periodic orbits can be found on this family of manifolds, shrinking to an 
equilibrium point at the HI point. 

A Hopf point of codimension 2 is a Hopf point where l\ vanishes. It is 
called transversal if rj — and l\ = have transversal intersections, where 
1] = i](fi) is the real part of the critical eigenvalues. In a neighborhood 
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of a transversal Hopf point of codimension 2 — H2 point, for concision — 
with I2 7^ the dynamic behavior of the system (fT4"l) . reduced to the family 
of parameter-dependent continuations of the center manifold, is orbitally 
topologically equivalent to 

w' — (77 + ito )w + rw\w\ 2 + l2iu\w | 4 , 

where 77 and r are unfolding parameters. See [Bj. The bifurcation diagrams 
for l 2 7^ can be found in [6], p. 313, and in [Tij . 

A i/op/ point of codimension 3 is a Hopf point of codimension 2 where 

1 2 vanishes. A Hopf point of codimension 3 point is called transversal if 
77 = 0, l\ = and Z2 = have transversal intersections. In a neighborhood 
of a transversal Hopf point of codimension 3 — H3 point, for concision — 
with Z3 7^ the dynamic behavior of the system ([ED , reduced to the family 
of parameter-dependent continuations of the center manifold, is orbitally 
topologically equivalent to 

w' = (77 + ito )w + rw\w\ 2 + uw\w\ A + lzw\w\ 6 , 

where rj, r and v are unfolding parameters. The bifurcation diagram for 

13 can be found in Takens [14] and in [13] . 

A Hopf point of codimension ^ is a Hopf point of codimension 3 where 
I3 vanishes. A Hopf point of codimension 4 is called transversal if 77 = 0, 
h — 0, li = and h = have transversal intersections. In a neighborhood 
of a transversal Hopf point of codimension 4 — H4 point, for concision — 
with Z 4 7^ the dynamic behavior of the system (Tl41) . reduced to the family 
of parameter-dependent continuations of the center manifold, is orbitally 
topologically equivalent to 

w ' = (v + ivo)w + tw\w\ 2 + vw\w\ 4 + aw\w\ 6 + Uw\w\ 8 , 

where 77, r, v and a are unfolding parameters. 
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Theorem 3.2 Suppose that the system 

x = /(x, fi), x = (x, y, 2), /i = (/3, a, ft, e) 
/ias £/ie equilibrium x = for fi = wi/i eigenvalues 

A 2 ,3(/i) = ± iw(^), 

where u(0) = ujq > 0. For /i = £/ie following conditions hold 

v(o) = o, h(0) = 0, / 2 (0) = 0, / 3 (0) = 0, 

where h(jJ>), and h(fi) are the first, second and third Lyapunov co- 

efficients, respectively. Assume that the following genericity conditions are 
satisfied 

1- £4(0) 7^ 0, where / 4 (0) is the fourth Lyapunov coefficient; 

2. the map \i — > (77 (>u) , h(n), hil^), hit 1 )) is regular at /1 = 0. 

Then, by the introduction of a complex variable, the above system reduced to 
the family of parameter- dependent continuations of the center manifold, is 
orbitally topologically equivalent to 

w' = (77 + iuo)w + tw\w\ 2 + uw\w\ A + aw\w\ 6 + l4w\w\ 8 

where r\, t, v and a are unfolding parameters. 



4 Hopf bifurcations in the WGSS 

The following theorem was proved by the authors in 



Theorem 4.1 Consider the four-parameter family of differential equations 
|3J). The first Lyapunov coefficient at the point (Q) for parameter values 
satisfying (TT3|) is given by 

1 (a \ Gi(/3,a,«) 

'' Wv) =-^( £!+s « + m' (35) 
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Figure 2: Signs of the first Lyapunov coefficient for system ([5]). 

where 

Gi(/3, a, k) = -3 + 5k/3 - (a 2 - 5)/? 2 + /t(a 2 - 7)/? 3 - 

2aV/3 4 - (a 4 - 2aV)/3 6 + a V 7 . (36) 

J/ G*i zs different from zero then the system (f3p has a transversal Hopf point 
at P for e = e c . More specifically, if (/3,a,K,) G S U U and e = e c then 
the system |3P has an HI point at P ; if (/3, a, k) G S and e = e c then 
the HI point at Pq is asymptotically stable and for each e < e C! but close 
to e c , there exists a stable periodic orbit near the unstable equilibrium point 
Poi if {fl, ct-i k) G U and e = e c then the HI point at P is unstable and for 
each e > e c , but close to e c , there exists an unstable periodic orbit near the 
asymptotically stable equilibrium point Pq. See Fig{^ 
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Figure 3: Signs of the second Lyapunov coefficient for system ©. 

Theorem 4.2 For the four-parameter family of differential equations |3]) 
there is unique point Q = a, k, e c ), with coordinates 

(3 = 0.93593 a= 1.02753 ...,k = 0.90164 e c = 0.73522 

where the surfaces l\ = 0, \<i = and I3 = on the critical hypersurface 
intersect and there do it transver sally. Moreover, the codimension 4 Hopf 
point at Pq is asymptotically stable since U(Q) < 0. More specifically, if 
(P, a, k) G Si U £2 U U\ and e = e c then the system |3P has an H2 point at 
Pq; if W, ct, k) G S% U ^2 and e = e c then the H2 point at Pq is asymptotically 
stable; if (/3, a, k) G U\ and e = e c then the H2 point at Pq is unstable. Along 
the curves C\ and Ci = C21 U C22 U {Q} of Fig\4\ h vanishes. If (/3, a, k) G 
C\ U C21 U C22 (see Fig and e = e c then the four-parameter family of 
differential equations ([3]) has a transversal Hopf point of codimension 3 at Pq; 



20 



if (J3, a, k) e G\ U C22 awe? £ = s c then the H3 point at Pq is asymptotically 
stable and the bifurcation diagram for a typical point H is draw in Fig 0|- 
if (/?, a, k) G C21 and e = e c then the H2 point at Pq is unstable and the 
bifurcation diagram for a typical point G can be found in JT%$ . 

Computer assisted Proof. The algebraic expression for the second Lya- 
punov coefficient can be obtained in [16]. This is too long to be put in print. 
The surface where the second Lyapunov coefficient vanishes is illustrated in 
Fig. El 



at 








Figure 4: Surfaces l\ = and I2 = and the intersection curves. 

The intersections of the surfaces li = and I2 = determine the curves 
C\ and Ci (see FigH]). The signs of the second Lyapunov coefficient on the 
surface l\ = complementary to the curves C\ and C 2 , that is on Si U S 2 U U\ 
(see Fig. [5]), are the following: Z 2 is negative on S\ U £2 and is positive on U\ 
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and they can be viewed as extensions of the signs of the second Lyapunov 
coefficient at points on the curve determined by the intersection of the surface 
li = and the plane k = studied by the authors in [12]. The bifurcation 
diagram for a typical point G where /3(G) > can be viewed in [12]. In Fig 
M and [7] are illustrated the bifurcation diagrams for a typical point H where 
h(H) < 0. 








Figure 5: Signs of l\, I2 and £3. 

The point Q is the intersection of the surfaces l\ = 0, I2 — and /3 = 0. 
The existence and uniqueness of Q with the above coordinates has been 
established numerically with the software MATHEMATICA 5. 
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Figure 6: Bifurcation diagram for a typical point H where h(H) < 0. 



For the point Q take five decimal round-off coordinates (3 = 0.93593, a = 
1.02753, k = 0.90164 and e c = 0.73522. For these values of the parameters 
one has 



p = H/2, 0.27041 - 0.54618*, 0.40395 + 0.20000z) , 
q = (-i, 0.36401, 0.99407), 
h n = (-2.65769,0,0.19650), 
h 20 = (-4.11029 - 0.18429i, 0.13416 - 2.9924H, 0.09159 - 3.36395i) , 
h 30 = (-3.63589 + 23.03616i, -25.15645 - 3.97054i, -18.16113 - 1.69167*) , 



h 2 i = (3.24775 + 1.67247z, -4.52694 + 1.18222*, 4.85950 + 3.71541z) , 



G 



3.91814*, 



(37) 



h. 




160.39204 + 51.10539* 
-74.41230 + 233.53975*' 
-25.03366 + 127.34049*' 



) 
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69.44664 - 38.56274i, 25.90851 - 2.24484i, 36.10391 - 65.85524i) , 
h 22 = (-64.50829, 0, 10.76131) , 

G 32 = -153.21726i, (38) 

/ 702.48693 - 1263.93346i \ 
h 50 =l 2300.44688+ 1278.5751 \% , 
\ 1054.20770 + 363.36145i / 

/ 178.24934 + 273.6678H \ 
h 32 =\ -233.17715 + 26.70966i , 
\ 395.89053 + 272.77265i / 

/ -521.71430 + 1074.2612H \ 
h 41 = -631.58388 - 484.25803i 

\ -865.10385 -413.20000i / 

/ -10130.73267- 9995.21750i \ 
h m = 21830.38995 - 22126.36639i , 
\ 5429.65950 - 9557.27148i / 

/ 14227.43860 + 8237.49829i \ 
h 51 = -9991.87299 + 14431.55078i , 
\ -5753.08267+ 11280.54380i / 

/ -4351.45992 - 4936.33553i \ 
h 42 = 2272.08822 + 1527.90723i 

\ 4841.97866 - 5445.36779i / 

h 33 = (-5969.63958, 0, 1764.47230) , 

/ -146941.54096 + 63522.80004i \ 
h 70 =\ -161862.28504- 374421.36634i , 
\ -86069.40319 - 83969.45215i / 

/ 140223.18890- 184094. 16057i \ 
h ei =i 260780.07852 + 213929.28545i , 
\ 151116.49070 + 92225.27059i / 

/ -105557.32750 + 127994.80577i \ 
h 52 =\ -41289.02476- 79039.91108i , 
\ -106857.14273 - 88122.45467i / 
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G 43 = -22328.21224i. (39) 

/ 26579.27090 + 62051. 16515z \ 
h 43 =l -36944.56779 + 2499.10743i , 
\ 78144.32459 + 54070.14624i / 

/ -247681.58290 + 2173895.03048i \ 
h m = -6330624.44741 - 721276.35507i , 
\ -1324248.15135 + 594661. 38331i / 

/ -2230744.30930- 2511854.85381* \ 
hn=\ 4663683.99275- 4038564.7541 li , 
\ 1618564.33911 - 2037646. 14488i / 

/ 2540059.79128 + 2277848.86298i \ 
h 62 =\ -2385453.21697+ 1869088.06376i , 
\ -1708253.47087 + 2025268.530342 / 

/ -633499.15640- 1125590.51413i \ 
hsa = 390598.08062 + 466226.40735i 

\ 1219484.73373 - 1101283.41903i / 

h 4A = (-1118100.12194,0.00138,546721.10946), 

G 54 = -22071.41115 - 5991090.521191 (40) 
From (EE}, (JHDD, ([32]), (El, (E7J), (EHD, (GEO) and gOl) one has 

h(Q) = 0, l 2 (Q) = 0, h(Q) = 0, U(Q) = -!- Re G 54 = -7.66368. 

zooU 

The calculations above have also been corroborated with 100 decimals 
round-off precision performed using the software MATHEMATICA 5 [17J. 
See pS]. 

Some values of (a, (3, k) G Ci U C 2 as well as the corresponding values of 
l s (a,/3,K) are listed in the tables below. The calculations leading to these 
values can be found in [16J. 
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Figure 7: Bifurcation diagram for a typical point H\. See Fig. [H 



K 


a 


P 


h(a, (3, k) on C\ 


0.45 


0.33319 


0.72216 


-0.91310 


0.5 


0.42968 


0.71770 


-0.92567 


0.55 


0.50934 


0.71257 


-0.88152 


0.6 


0.57913 


0.70665 


-0.82064 


0.65 


0.64241 


0.69983 


-0.75810 


0.7 


0.70113 


0.69201 


-0.70006 


0.75 


0.75659 


0.68309 


-0.64900 


0.8 


0.80972 


0.67302 


-0.60580 


0.85 


0.86120 


0.66177 


-0.57054 


0.9 


0.91154 


0.64940 


-0.54288 


0.95 


0.96114 


0.63600 


-0.52217 
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K 




B 




n 


85050 


86898 








87760 


D 46294 


0.3 


0.93123 


0.88397 


0.50684 


0.4 


0.95511 


0.89159 


0.55538 


0.5 


0.97602 


0.90042 


0.60637 


0.6 


0.99330 


0.91029 


0.65253 


0.7 


1.00674 


0.92071 


0.66963 


0.8 


1.01697 


0.93045 


0.56860 


0.9 


1.02731 


0.93592 


0.01665 


0.92 


1.03020 


0.93585 


-0.20674 


0.98 


1.04319 


0.93201 


-1.09289 



The gradients of the functions h, h and I3, given in (1281) . (|30l) . ( 1321) at 
the point Q are, respectively 

(-0.46264, 0.13437, -0.97565), (-12.44701, 2.66791, -19.19345), 

(-266.77145, 41.80505, -372.84969). 

The transversality condition at Q is equivalent to the non-vanishing of the 
determinant of the matrix whose columns are the above gradient vectors, 
which is evaluated gives —33.31133. 

The main steps of the calculations that provide the numerical evidence 
for this theorem have been posted in [16J. 



5 Concluding comments 

This paper starts reviewing the stability analysis which accounts for the 
characterization, in the space of parameters, of the structural as well as 
Lyapunov stability of the equilibrium of the Watt Governor System with a 
Spring, WGSS. It continues with recounting the extension of the analysis to 
the first order, codimension one stable points, happening on the complement 
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of a surface in the critical hypersurface where the eigenvalue criterium of Lya- 
punov holds, as studied by the authors [15] . based on the calculation of the 
first Lyapunov coefficient. Here the bifurcation analysis at the equilibrium 
point of the WGSS is pushed forward to the calculation of the second, third 
and fourth Lyapunov coefficients which make possible the determination of 
the Lyapunov as well as higher order structural stability at the equilibrium 
point. See also © Ej , [I] and [2] . 

The calculations of these coefficients, being extensive, rely on Computer 
Algebra and Numerical evaluations carried out with the software MATHE- 
MATICA 5 [17]. In the site [16] have been posted the main steps of the 
calculations in the form of notebooks for MATHEMATICA 5. 

With the analytic and numeric data provided in the analysis performed 
here, the bifurcation diagrams are established along the points of the surface 
where the first Lyapunov coefficient vanishes. Pictures [6] and [7J provide a 
qualitative synthesis of the dynamical conclusions achieved here at the pa- 
rameter values where the WGSS achieves most complex equilibrium point. 
A reformulation of these conclusions follow: 

There is a "solid tongue" where three stable regimes coexist: one is an equilib- 
rium and the other two are small amplitude periodic orbits, i.e., oscillations. 

For parameters inside the "tongue" , this conclusion suggests, a hysteresis 
explanation for the phenomenon of "hunting" observed in the performance 
of WGSS in an early stage of the research on its stability conditions. Which 
attractor represents the actual state of the system will depend on the path 
along which the parameters evolve to reach their actual values of the param- 
eters under consideration. See Denny [3] for historical comments, where he 
refers to the term "hunting" to mean an oscillation around an equilibrium 
going near but not reaching it. 

Finally, we would like to stress that although this work ultimately focuses 
the specific three dimensional, four parameter system of differential equations 
given by (jSJ), the method of analysis and calculations explained in Section 
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[3] can be adapted to the study of other systems with three or more phase 
variables and depending on four or more parameters. 
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